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We introduce a quantum dimer model on the hexagonal lattice that, in addition to the standard 
three-dimer kinetic and potential terms, includes a competing potential part counting dimer-free 
hexagons. The zero-temperature phase diagram is studied by means of quantum Monte Carlo simu¬ 
lations, supplemented by variational arguments. It reveals some new crystalline phases and a cascade 
of transitions with rapidly changing flux (tilt in the height language). We analyze perturbatively 
the vicinity of the Rokhsar-Kivelson point, showing that this model has the microscopic ingredients 
needed for the “devil’s staircase” scenario [E. Fradkin et al. Phys. Rev. B 69, 224415 (2004)], and 
is therefore expected to produce fractal variations of the ground-state flux. 

PACS numbers: 05.30.-d, 05.30.Rt, 05.50.-|-q, 75.10.Jm, 


The study of hard-core dimer coverings has a long his¬ 
tory. From the mapping to Pfafiians and determinants 
by Kasteleyn [1, 2], the solution of two-dimensional Ising 
models [3], the height representation and its continuum 
limit [4], or the connection to the Coulomb gas and con¬ 
formal field theory [-5, 6], dimer models have found nu¬ 
merous applications in various fields of statistical physics. 
Motivated by the physics of resonating valence bond sys¬ 
tems, Rokhsar and Kivelson (RK) [7] added quantum dy¬ 
namics to the dimer model, leading to the so-called quan¬ 
tum dimer model (QDM), which later led to tractable 
models with rich phase diagrams closely related to lat¬ 
tice gauge theories [8]. Importantly, QDMs appeared in 
different contexts when describing the dynamics in a con¬ 
strained low-energy manifold, such as in frustrated Ising 
models in weak transverse fields [9]. QDMs also gained 
a new dimension with the discovery of liquid phases 
with topological order in nonbipartite lattices [10, 11], 
where they shed some light on the long-sought resonat¬ 
ing valence bond liquids. This field also benefited from 
recent progress in making quantitative connections be¬ 
tween spin-1/2 Heisenberg magnets with quantum disor¬ 
dered ground states and QDMs [12, 13]. 

In most QDMs studied so far, a kinetic term (associ¬ 
ated with on-plaquette dimer flips) competes with a diag¬ 
onal term proportional to the number of such “flippable” 
plaquettes. When the kinetic and the potential terms 
are equal at the so-called RK point, the ground states 
are exactly known [7]. In the height language, appropri¬ 
ate for bipartite lattices, such a RK point corresponds 
to a transition from a “flat” phase to a maximal slope 
phase [14]. A richer behavior is however expected near 
that point for more generic interactions between dimers 
[15, 16]. In particular, within a field theoretic approach, 
a devil’s staircase of commensurate and incommensurate 


phases is predicted [15-17], corresponding to a fractal tilt 
variation as a function of the Hamiltonian parameters. 

In this Letter, we show that a natural generalization 
of the hexagonal lattice QDM [18, 19] provides a micro¬ 
scopic model with this phase structure. We analyze the 
two-parameter phase diagram spanned by the standard 
potential term counting flippable plaquettes and another 
term counting dimer-free plaquettes. The model is stud¬ 
ied perturbatively near the RK point and with quan¬ 
tum Monte Carlo (QMC) simulations elsewhere, sup¬ 
plemented by variational arguments. We observe a se¬ 
quence of closely spaced phase transitions with a gradual 
change of the flux density and crystalline structures with 
strongly varying unit cell sizes in agreement with the sce¬ 
nario of Refs. [15, 16]. 

Model. —Let us consider a QDM with the standard ki¬ 
netic term and four potential terms: 
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H = - ([/_» (/->[ -b H.c.) -b ^ Vjfij, (1) 

h j=0 

where the operator nj counts the total number of hexag¬ 
onal plaquettes with j dimers (called a j-plaquette). Be¬ 
cause of the two sum rules [19, 20] no + ni+ 722+713 = N 
and 2ho + fii — 223 = 0 , these potential terms are not 
independent and we hence choose to keep only fiQ and 
723 . Also, we denote densities pj = {nj)/N in the form 
p = (po, pi, p 2 , ps) and fix t = 1, unless specified differ¬ 
ently. The model studied by Moessner et al. [18] has 
vq = 0, while the two models (vq = ± 1 ,U 3 = 0) are rele¬ 
vant for Ising string nets [21]. We study rectangular clus¬ 
ters with periodic boundary conditions and N = x Ly 
hexagonal plaquettes. 

Our analysis relies on the notion of flux: dimer cover¬ 
ings can be grouped into topological sectors [20] labeled 
by two integer fluxes (F/, F/), which are invariant under 
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FIG. 1. Schematic phase diagram from QMC simulations 
{Lx = Ly = 60). The (vojVs) plane is divided into five re¬ 
gions: a staggered phase with the maximal flux (/ = 2), 
the star and the plaquette phases (/ = 0), the S 2 phase 
(/ = 1/2), and the fan region, containing a cascade of flux 
sectors 1/2 < / < 2. The plaquette color indicates the dimer 
density (same scale as Figs. 2 and 3). 

local dimer moves. As discussed below, for ground states, 
one of the two fluxes is zero and we can restrict ourselves 
to Fx = 0 and work with / •= Fy/Ly > 0. 

Classical limit. —Let us consider the classical limit 
t = 0. Setting uo = sina, U 3 = cos a, and defining 
ai = arctan(—2), ot2 = 7r/2 — ai, one finds three crys¬ 
tals as a is varied: (i) for a S [7r/2, cti], the three¬ 
fold degenerate staggered crystals (nonflippable config¬ 
urations) with maximum flux / = 2 , vanishing energy, 
and p = ( 0 , 0 , 1 , 0 ), (ii) for a S [ 01 , 0 : 2 ], the (threefold 
degenerate) star crystal in the / = 0 sector (Fig. 1) with 
p = (1/3,0,0, 2/3), (iii) for a G [ 02 ,7’'/2], a 12-fold degen¬ 
erate crystal [ 22 ] denoted S' 2 , within the / = 1/2 sector, 
with p = (0,1/2, 0,1/2). The point a = 7r/2 is highly 
degenerate, since any configuration without 0 -plaquettes 
is a ground state, and such states exist in all flux sec¬ 
tors. This degeneracy is lifted when t yf 0, leading to a 
nontrivial ground-state flux variation as discussed below. 

Phase diagram. —We studied the phase diagram with 
QMC simulations using the mapping to an Ising-type 
model described in Refs. [18-20]. Specifically, results 
displayed in Fig. 1 have been obtained for a torus with 
60 X 60 plaquettes, flux sectors / = 0, ^ 2 , in¬ 

verse temperature P = 9.6, and imaginary-time step 



FIG. 2. Left: a configuration of three strings and the cor¬ 
responding dimer covering, with 0- and 3-plaquettes. Left 
bottom and right: three variational classes of dynamically 
constrained strings, called S-, H-, and F-strings. F-strings 
are in a static zigzag configuration, iF-strings (F-strings) are 
allowed to fluctuate by one row in every second column (in ev¬ 
ery column). Arrows indicate the fluctuations of the strings, 
each corresponding to a 3-plaquette flip. Dimer densities are 
indicated according to the color scale of Fig. 3. For H- and 
F-strings, the shown dimer densities correspond to a super¬ 
position of the allowed configurations. 

AP = 0 . 01 . 

• / = 2. In this region, ground states are isolated stag¬ 
gered configurations with vanishing energy. The Hamil¬ 
tonian is positive definite in the upper right quadrant, 
and the / = 2 region also extends to a large part of the 
lower right quadrant, down to the boundary with the 
/ = 0 sector. 

• / = 0. The star and plaquette crystals found in this 
region also exist in the ua-only model [18, 19] and are 
separated by a first-order transition (dashed line). The 
star phase is adiabatically connected to the (threefold 
degenerate) crystalline configurations found for t = 0. 
The latter simultaneously maximize the number of 3- 
and 0 -plaquettes, and the star phase thus fills a large 
part of the {v^ < 0 ,vo < 0 )-quadrant and also extends 
into the neighboring quadrants. On the uq = 0 line, the 
star phase gives way to the plaquette phase through a 
first-order transition at V 3 = —0.228(2) [18, 19]. The 
plaquette phase is defined by continuity with the “ideal” 
plaquette state, which is an uncorrelated product of res¬ 
onating 3-plaquettes -I- 1\~/). In the vicinity of the 
RK point, as is already the case for H{t, vq = 0, V 3 ) [19], 
the large (diverging) correlation length makes it difficult 
to discriminate numerically between the star and pla¬ 
quette phases, hence the question mark in Fig. 1. This 
phenomenon is likely to be related to the U{1) regime 
observed in the square lattice QDM [23]. 

• f — tj^l. In most of this region, the system forms 
a 12 -fold degenerate crystalline phase, adiabatically con¬ 
nected to the S '2 configuration. 

• l/2</<2. This is the most interesting part of the 
phase diagram, which we call the fan region. To under¬ 
stand the flux variations taking place there, we recall that 
any dimer configuration can be represented equivalently 
as a configuration of nonintersecting strings on the hexag- 
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onal lattice [20]. For = 0, these are Ng = {2Ly — Fy)/3 
closed loops along the toroidal x direction of the lattice. 
Starting from the staggered dimer covering (/ = 2) dis¬ 
played in Fig. 1, on each string path, empty and covered 
edges alternate. The corresponding dimer covering is ob¬ 
tained by doing so-called loop updates, i.e., exchanging 
empty and covered edges along the string paths. Each 
string reduces the flux Fy by three units. In reverse, 
starting from an arbitrary configuration, the strings cor¬ 
respond to paths where dimer-free horizontal edges alter¬ 
nate with dimers on tilted edges (see Fig. 2). The number 
of 3-plaquettes along a string is maximized if it runs par¬ 
allel to one of the three edge orientations of the lattice. 
This is why, for U 3 < 1, strings are on average parallel to 
one of the edge orientations and why ground states are 
found in sectors with one vanishing flux quantum num¬ 
ber {Fx vanishes for strings winding in the x direction 
only). Strings can reduce their kinetic energy by oscillat¬ 
ing in the perpendicular direction, limited by the string 
noncrossing condition and by avoidance of 0 -plaquettes 
for large uq (see Fig. 2). 

When U 3 is decreased below 1, the staggered configu¬ 
ration is destabilized by string insertion. At low string 
densities (/ slightly below 2 ) strings are far apart and 
strongly delocalized. A reduction of U 3 causes an increase 
of p 3 , which is realized through a higher string density 
(decrease of the flux) and “stiffer” strings (reduced lat¬ 
eral motion). Each time a new string is added upon de¬ 
creasing U 3 , the increased compensates the energy cost 
associated with the higher degree of localization. When 
increasing uq for a fixed U 3 < 1 , configurations with more 
0 -plaquettes become less favorable such that string de- 
localization gets more restricted. At certain transition 
points, it becomes favorable to remove a string (flux in¬ 
crease), freeing some space for other strings to fluctuate 
more freely. When po becomes negligible, a further in¬ 
crease of Vq has no effect. This regime, where the isoflux 
lines become parallel, is equivalent to perturbing the (de¬ 
generate) classical point {t, vq, V 3 ) = ( 0 , 1 , 0 ) with a weak 
t and V 3 , where a fan like phase diagram similar to that 
described in Ref. [24] is expected. 

For / ^ 1 the average interstring distance is suffi¬ 
ciently low that the ground states are dominated by 
straight-string configurations. For generic fluxes, one 
expects complex correlated string states (some are de¬ 
scribed in Ref. [20]), but simple spatial structures in¬ 
volving horizontal chains of hexagons with higher den¬ 
sities of 3-plaquettes are also observed in some low-flux 
parts of the fan (see Fig. 3). These can be qualitatively 
understood in terms of the following typical configura¬ 
tions of strings that are dynamically constrained by the 
presence of neighboring strings: “5'-strings” are static 
zigzag configurations (corresponding to zigzag arrange¬ 
ments of 3-plaquettes, energetically favored at large neg¬ 
ative V3). With respect to such a reference configura¬ 
tion, “i7-strings” can fluctuate in every second column 
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FIG. 3. Dimer density per hexagon in the fan region, accord¬ 
ing to the color scale on the right. 


of hexagons, up and down by one row. “E-strings” are 
the most mobile among the three classes, and are allowed 
to fluctuate up and down by one row in every column as 
indicated by arrows in Fig. 2. At / = 0.8 and 1, for 
instance, we recognize periodic arrays oi FI- {F-) strings 
at distance d = 2.5 {d = 3) [25], as shown in Fig. 3. 
Importantly, no 0-plaquettes are generated if the above 
strings have minimum interstring distances of = 2 , 
= 2.5, df^p = 3, and df% = 2.75. These build¬ 
ing blocks are therefore appropriate to describe qualita¬ 
tively the large-uo and / ;< 1 part of the fan [ 20 ]. 

Finally, simple variational arguments provide approxi¬ 
mate expressions for the flux transition lines. For exam¬ 
ple, one can compute the energy change associated with 
the insertion of an i7-string in a perfect S 2 crystal {S- 
strings at distance 2 ), which corresponds to an infinites¬ 
imal increase of the flux density (due to the different 
^mm, gyg E-strings should be replaced by four i7-strings 
to keep the total system size constant [20]). This yields 
U 3 = — I for the transition towards the fan region at large 
Vg, in reasonable agreement with the numerics. 

As the interplay between 1:3 and vg is especially com¬ 
plex for low vg (when pg is not negligible), we analyzed 
the U 3 = 0 line with finer flux steps. Starting from very 
large vg the flux decreases (staying close to / = 0 . 8 ) down 
to Wo — 2.4 where it drops to / = 0. This flux drop is 
a generic feature of the interface with the / = 0 region. 
Toward the RK point the ground-state flux sectors get 
pinched, a feature that we now discuss. 

Perturbative analysis. —At the RK point, the ground 
states of all flux sectors are degenerate, and are equal- 
amplitude superpositions of all dimer configurations 
in the corresponding sector. At first order in vg/t 
and {v 3 — l)/t, the energy density in sector / reads 
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flux f angle e = atan[vo/{v3-1 )] 


FIG. 4. Perturbation theory near the RK point. Left: j- 
plaquette densities pj as functions of the flux density /. Right: 
ground-state flux density / as a function of the angle 6 that 
parametrizes the perturbation. The nontrivial region lies be¬ 
tween 7r/2 (transition out of the staggered phase) and 9i ~ 
1.84695, where / drops discontinuously from /i ~ 0.195654 to 
0. The transition from / = 0 to / = 2 occurs at 62 — 4.8268 
(not shown). 

e(/) ='yoPo(/) + - 1 )P 3 (/)- We compute the j- 

plaquette densities Pj{f) as expectation values of the op¬ 
erators fij (diagonal in the dimer basis) with respect to 
the unperturbed RK states, using an analytical transfer- 
matrix approach [20, 26]. Setting vq = sin 6* and ^3 — 1 = 
COS0, we minimize e(/) for each value of 9 to obtain f{9) 
as displayed in Fig. 4. A continuous variation of / is 
found in the interval 0 S [ 7 r/ 2 , 0i ~ 1.84695], which cor¬ 
responds to the fan region in the phase diagram of Fig. 1. 
Interestingly, / jumps discontinuously to zero at 9i. For 
0 6 [ 01,02 — 4.8268], the ground state is in the / = 0 flux 
sector, and it jumps to / = 2 for 0 G [ 02 , 7 r/ 2 ]. Note that, 
at this order, wave functions remain RK states, which are 
translation-invariant dimer liquids with algebraic corre¬ 
lations (for / < 2 ). 

Field theory. —To connect our perturbative and numer¬ 
ical results concerning the flux variations, let us turn to 
the height representation [4, 17, 27-29]. Dimer coverings 
are mapped to membranes embedded in a cubic lattice, 
whose average tilt is directly related to the flux [30]. In 
this language the QDM becomes a quantum roughening 
problem [17]. Long-distance properties are captured by 
taking the continuum limit of the height model and, in 
our case, the RK point is described by a massless Gaus¬ 
sian field theory [29]. Fradkin et al. [15] and Vishwanath 
et al. [16] discussed how the action is modified in the 
presence of generic perturbations, through a renormal¬ 
ization group (RG) analysis [31] predicting nonvanishing 
flux phases. A cubic interaction for the height, with three 
spatial derivatives, is the leading term favoring / ^ 0 . 
In our problem we observe that Vg induces a flux den¬ 
sity perpendicular to some edges of the hexagonal lattice. 
This implies that the sign of the corresponding coupling 


is negative in the notation of Ref. [15]. At this stage, the 
system would be gapless with a linear dispersion at small 
momenta. However, the site positions and the micro¬ 
scopic heights are both discrete and form a 3D lattice £. 
For the (coarse-grained) height field, potential terms that 
respect the symmetries of £ will be generated upon in¬ 
tegration over the short-distance fluctuations. They can 
be written as V(h,rJ = EK=(iyo,£)G£* 
where the sum runs over the reciprocal lattice vectors of 
£. When the average flux (tilt) is commensurate with the 
lattice, it corresponds to some reciprocal lattice vector 
K and the associated locking term Vk is then asymptot¬ 
ically relevant in the RG [15], leading to gapped crystals. 
However, as explained in Ref. [15], these gaps can be¬ 
come exponentially small in 1// close to the RK point. 
Since crystals for rational fluxes with small denomina¬ 
tors are more stable, their range of attraction in the RG 
is larger compared to others and, for the phase diagram 
close to the RK point, one thus expects a fractal suc¬ 
cession of commensurate phases — a “devil’s staircase”. 
At the smaller fluxes, stronger quantum fluctuations can 
outweigh locking terms and impose irrational flux den¬ 
sities such that gapless incommensurate structures are 
possible. 

Conclusion .—The extended QDM (1) is the first can¬ 
didate for a microscopic realization of the “Gantor decon¬ 
finement” scenario, which predicts that a fractal succes¬ 
sion of flux sectors occurs near the RK point. Whether 
the flux varies continuously, in a fractal way, or assumes 
only a finite number of values [32] is impossible to an¬ 
swer with QMG simulations. Indeed, although we can 
simulate large lattices, available flux sectors correspond 
to a small set of rational values. Additionally, intrasec¬ 
tor gaps become very small near the RK point and render 
simulations difficult. However, the fact that all flux sec¬ 
tors for 1/2 < / < 2 occur in the QMG results and the 
width variations of the corresponding regions in the phase 
diagram plead in favor of the realization of a fractal in 
the thermodynamic limit. 

Finally, let us note that flux sequences found here can¬ 
not occur for square lattice models with single-plaquette 
Hamiltonians. In that case, the sum rule no = n 2 makes 
any QDM with potential terms Vjhj equivalent to the 
original RK model, which lacks intermediate-flux phases. 
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Supplemental Material for “Phase 
Diagram of an Extended Quantum Dimer 
Model on the Hexagonal Lattice” 


Appendix A: Topological flux quantum numbers 

The hexagonal lattice is bipartite: all the neighbors of 
a site in the “even” sublattice belong to the “odd” sub¬ 
lattice, and vice versa. Let us recall that, with periodic 
boundary conditions, the set of all dimer coverings on 
a bipartite lattice breaks into topological sectors which 
are stable under any local dimer move (including single- 
plaquette flips). These sectors can be labeled by a pair of 
flux quantum numbers {Fj.,Fy). Locally, a “magnetic” 
field B can be defined as follows. As shown in Fig. 5, an 
empty edge carries one unit of magnetic field, oriented 
from the even to the odd sublattice. And each edge oc¬ 
cupied by a dimer carries two field units going from the 
odd to the even sublattice. We shall count the magnetic 
flux going through oriented closed paths, the smallest 
one being a small circle surrounding a single site of the 
hexagonal lattice. Due to the constraint that every site is 
reached by exactly one dimer, the flux going through such 
a small circle vanishes - the lattice divergence of the mag¬ 
netic held is zero at each site {divB = 0). Correspond¬ 
ingly, the hux through any contractible loop, such as C 
in Fig. 5, also vanishes. Indeed, since divB = 0, one can 
deform the loop until it reduces to the smallest one. On 
the other hand, the hux through non-contractible loops 
can be non-zero. Again, any local dimer rearrangement, 
including the hip term of the Hamiltonian, conserves the 
hux through such loops. Since, for the torus, there exist 
two independent non-contractible loops, this leads to a 
pair of huxes F^ and Fy characterizing hip-disconnected 
topological sectors. The star and plaquette crystals, dis¬ 
cussed in the main text, are found in the zero-hux sector, 
while the staggered dimer conhgurations are found in the 
maximal hux sectors, each consisting of a single (isolated) 
conhguration. 


Appendix B: String representation and fluxes 

Dimer coverings can equivalently be represented as 
conhgurations of nonintersecting strings that form closed 
loops on the lattice [19, 33]. Specihcally, one can start 
from the staggered dimer conhguration depicted on the 
left of Fig. 5 (horizontal dimers only) and choose a certain 
conhguration of closed nonintersecting strings such that, 
on the path of every string, empty and occupied edges 
alternate. Now, the corresponding dimer conhguration 
is obtained by doing so-called loop updates along the 
string paths. These consist in exchanging, on the chosen 


paths, empty and occupied edges. The winding num¬ 
bers of these strings are in direct correspondence with 
the hux quantum numbers {Fx,Fy) introduced above. If 
we consider, for example, the case with F^ = 0 that is 
the relevant sector for the main part of the paper, every 
dimer covering with hux Fy corresponds to a conhgura¬ 
tion of Ns strings that encircle the torus in x direction 
such that 

Fy = 2Ly-3Ns ^ f = Fy/Ly=2-ms/Ly. ( 2 ) 


Appendix C: Height representation 

Dimer coverings of a bipartite lattice can equivalently 
be represented using an integer height associated with 
each plaquette (dual lattice sites) [4, 34]. On the hexag¬ 
onal lattice, turning clockwise around a site of the even 
sublattice, the height h{x) changes by -|-1 when crossing 
an empty edge, and by -2 when crossing a dimer (respec¬ 
tively -1 and +2 around a site of the odd sublattice), 
as depicted in Fig. 6. The magnetic held B, dehned 
above, is perpendicular to the slope in this height rep¬ 
resentation. So, conhgurations conhgurations with zero 
hux have a vanishing average slope in the height language 
(“hat” conhgurations). Similarly, conhgurations with a 
large hux correspond to a large slope. This maps the 
dimer covering problems onto faceting problems for sur¬ 
faces in three dimensions [17]. It is indeed conventional 
to represent dimer coverings on the hexagonal lattice as 
“stacks of cubes” (see for instance Fig. lb in [17]). The 
height h introduced above then corresponds to the po¬ 
sition of the surface of the cubes, after projection onto 
the (1,1,1) axis of the underlying cubic lattice. In this 
language, staggered dimer regions appear as smooth (but 
tilted) surfaces, normal to the (1,0,0), (0,1,0) or (0,0,1) 
direction. On the other hand, (static) iS-strings at dis¬ 
tance d from each other correspond to steps separating 
(tilted) terraces of width d. The star crystal corresponds 
to a microscopically corrugated surface with a vanishing 
average slope. 


Appendix D: Dimer sum rules 

Dimer coverings on a tiling by dehnition satisfy a sim¬ 
ple constraint - each vertex is reached by exactly one 
dimer. Hence, dimer coverings are constrained by simple 
sum rules, associated with Euler-Poincare and Gauss- 
Bonnet relations for tilings on compact surfaces [35]. 
For a given covering, let Nd denote the total number 
of dimers, nj the number of plaquettes covered with j 
dimers, and N and V the total number of plaquettes and 
vertices, respectively. Calling jmax the maximum num¬ 
ber of dimers that can sit on a plaquette, this gives the 


7 



FIG. 5. (Color online) Left: The perfectly staggered state (/ = 2) and a configuration of two strings on a torns with 8x5 
plaquettes. Right: The dimer covering that corresponds to the string configuration on the left is obtained by doing loop updates 
on the strings’ paths, exchanging empty and occupied edges. Dimers carry two units of flux and empty edges one unit of flux 
in the indicated directions. Due to the dimer constraints, the flux through contractible loops such as C is always zero and the 
two nontrivial flux quantum numbers and Fy for the 2D torus are invariant under local rearrangements. Only global loop 
updates, such as those indicated, change F^ and Fy. The color of each plaquette indicates its dimer density, according to the 
scale shown on the right (same as in Fig. 3 in the main text). 


(a) (b) (c) 



FIG. 6. Height mapping for the hexagonal dimer problem. 
An integer height is associated with each plaquette according 
to the rule given in the text, (a) A staggered dimer config¬ 
uration corresponds to a maximal slope, and therefore the 
maximal flux sector, (b) Strings (in this case an S-string) 
reduce the slope, (c) In the height representation, the star 
crystal has (on average) zero slope. 


two sum rules 


Appendix E: Classical phase diagram 

Let us compute the classical ground states [t = 0) 
for arbitrary values of vg and vg. To determine the j- 
plaquette densities p = {po, pi, p 2 , ps) that minimize the 
energy, we first introduce a parametrization of the acces¬ 
sible phase space. According to the sum rules 

3 

2/Oo + Pi - P3 = 0 and ^ pj = 1, (5) 

j=o 

the physical states are restricted to a triangular region in 
the (po,Pi,P 3 )-space, formed by the origin O = (0,0,0) 
and the points A = (1/3,0, 2/3) and B = (0,1/2,1/2). 
We parametrize a generic point P inside that triangle by 


(?max Jmax 

j • rij = 2Nd = V and rij = N. (3) 
j=o j=o 

The first follows from the local dimer constraint which 
implies that V = 2Nii. The second rule simply expresses 
that the total number of plaquettes is N. For a regular 
tiling, i.e., a tiling with constant coordination number c, 
the numbers of plaquettes and vertices obey N = V (c/2— 
1) on a torus. 

For the hexagonal lattice, we have c = 3, jmax = 3, 
and V = 2N, leading to 

3 

2nQ + m — rig = 0 and rij = N. (4) 

3=0 

Note that, on average, plaquettes carry two dimers. For 
the square lattice, c = 4, j^ax = 2, and V = N, lead¬ 
ing to the 712 = n-o result noted in the conclusion. Note 
finally that tilings with fixed boundaries can also be ana¬ 
lyzed along the same line, by properly entering additional 
boundary terms. 


, s 1 — s 1 s 

P= I + “ 


with (r,s)G [0,1]. (6) 


Now, with vg = sin a and vg = cos a, the energy per 
plaquette E(a) reads 


E{a) 


PoCo + P3C3, 

s ( . cosaN 

-(s.„a+—) 


+ - cos a 


( 7 ) 

( 8 ) 


We now seek for the minimum of E{pt) in terms of r and 
s. Clearly, the sign of (sin a-|-) decides whether 
s = 0 or 1 for the ground state configuration. With 
ai = arctan(—2) ~ —63.4°, a 2 = 7r/2 —ai and clockwise 
rotation, one obtains the three regions given in the main 
text: 

(i) a € [7r/2,ai]: E(a) > 0 and the ground state en¬ 
ergy is minimized (and vanishes) throughout this interval 
when 7- = 0, leading to p= (0, 0,1, 0). This corresponds 
to the staggered states (nonflippable configurations) in 
the maximal / = 2 flux sector. When a = 7r/2, any con¬ 
figuration satisfying s = 0, hence the full segment OB, 
also defines a ground state. Such configurations can be 











found in every flux sector. At the opposite end of this 
angular sector, when a = ai, all configurations falling in 
the segment OA, therefore s = 1 and r € [ 0 , 1 ], also have 
a vanishing energy. Again, such configurations exist in 
every flux sector. 

(ii) a S [ai, 02 ]: The threefold degenerate ground state 
is the star crystal, corresponding to the point A with 
p = (1/3,0,0, 2/3). It belongs to the / = 0 sector (Fig. 1, 
main text). When a = 02 , all configurations falling in 
the segment AS, with r = 1 and s G [ 0 , 1 ], have mini¬ 
mal energy. Such configurations can be found at least in 
sectors / G [0,1/2]. 

(hi) a G [a 2 , 7 r/ 2 ]: The ground state is a 12- 
fold degenerate crystalline state, denoted by S' 2 , with 
p = (0,1/2,0,1/2). It belongs to the f = lj2 sector. 

For this classical limit, the locations of ground states in 
the (poj Pi, P 3 )-space show a nice “dual” relation with re¬ 
spect to the circle (angle a) that parametrizes the Hamil¬ 
tonian. In the (uq, U 3 )-plane, the classical phase dia¬ 
gram has three angular sectors separated by the values 
( 7 r/ 2 , ai, 02 ). The set of classical configurations defines 
a convex region bounded by the triangle (O, A, B) in the 
(po, Pi,P 3 )-space. As described above, the ground states 
lie on the boundary of that triangle. Therefore, the con¬ 
tinuous a intervals map onto the triangle’s vertices, while 
the three singular values of a are mapped onto the tri¬ 
angle’s edges. 

Appendix F: Details about Monte Carlo simulations 

The numerical method used in this work has been de¬ 
tailed in Ref. [19], and we therefore only briefly summa¬ 
rize it here. As done by Moessner, Sondhi, and Chandra 
[18], the 2D quantum dimer model on a hexagonal lattice 
can be studied by first mapping it to an antiferromag¬ 
netic 2D quantum Ising model on the (dual) triangular 
lattice, comprising diagonal six-spin interactions. The re¬ 
sulting model can be studied efficiently using world-line 
quantum Monte Carlo [36] by approximating its partition 
function and observables by those of a classical 3D Ising- 
type model on a stack of triangular 2D lattices (quantum- 
classical mapping). We speed up the Monte Carlo simu¬ 
lation of the classical 3D model through suitable cluster 
updates. 

The equivalence between dimer and spin models is a 
delicate issue for two reasons. First, as we are free to 
choose the orientation of some reference spin, a given 
dimer configuration corresponds to two spin configura¬ 
tions that differ by a global spin flip. Dimer configu¬ 
rations therefore correspond to the spin-flip symmetric 
sector of the spin model. One can nevertheless simu¬ 
late the full spin model in the Monte Carlo as, based 
on the Perron-Frobenius theorem, it can be shown that 
the global ground state is always in the spin-flip sym¬ 
metric sector. Second, half of the topological sectors of 



FIG. 7. Points in the wo — wa plane where the simulations 
were carried out. The corresponding ground state fluxes (see 
color scale) were used to determine the flux transition lines 
shown here (grey lines) and in Fig. 1 (main text). 

the dimer configurations correspond to periodic bound¬ 
ary conditions for the spins, and half of the sectors corre¬ 
spond to anti-periodic boundary conditions in the Ising 
model. 

In the present work, we have only used periodic bound¬ 
ary conditions for the Ising spins, and systematically 
compared the available flux sectors in one direction (say 
Fy), while the other {Fx) is kept zero, by explicitly con¬ 
structing an appropriate initial spin configuration. In 
doing so, we have assumed that states in each flux sec¬ 
tor are flip connected. It is indeed generally believed, 
that the local dynamics are ergodic in each topological 
sector, besides those of maximum flux (see Prop. 2.3 in 
Ref. [37]). Most of the simulations were carried out for 
a rectangular cluster with 60 x 60 plaquettes, as noted 
in the main text, allowing for 20 different flux densities 
/y, equally spaced by steps of 5/ = 0.1. The points in 
the vq — V 3 plane that were investigated to determine the 
flux transition lines in Fig. 1 (main text) are displayed 
in Fig. 7. Other system sizes, ranging from 56 x 55 to 
120 X 120 plaquettes, were also used, mainly to study in 
more detail the transitions along the U3 = 0 line (cor¬ 
responding to the cluster of points for 2.5 < uq < 4 in 
Fig. 7). We verified in particular that the energy density 
is the same for two systems with different sizes but same 
flux density. 

For the QDM at uq = 0, the transition between the 
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star and plaquette phase occurs in the zero-flux sector 
and can be detected using magnetization variances of the 
corresponding spin model [18, 19]. In the more generic 
case with nonzero Uq, studied in this work, many tran¬ 
sitions between different flux sectors occur. These must 
therefore be detected by comparing the energies of dif¬ 
ferent sectors. As detailed in Ref. [19], the energy can 
be evaluated using imaginary-time spin-spin correlators 
We scanned the {vq,v^) plane in order to de¬ 
termine the ground state flux. In addition, the expecta¬ 
tion values of several spin and dimer observables (dimer 
densities and correlations) were measured to character¬ 
ize each phase. In particular, plots of the average dimer 
occupancies (for each plaquette) are used to visualize the 
spatial/crystalline structures and compare stripe-like or¬ 
ganizations in the “fan” region. 

Appendix G: Perturbation aronnd the RK point 
Classical transfer matrix and free fermions 

In this section we compute the dimer density expec¬ 
tation values pq and pa in the RK ground state, as a 
function of the flux density /. RK states being equal 
amplitude superpositions of all covering in a given flux 
sector, po and pa are also the densities of a classical sta¬ 
tistical dimer problem at infinite temperature. We solve 
the later using a transfer matrix method. 


y 



FIG. 8. Brickwall representation of the hexagonal lattice. 
The fermions (blue circles) live on the vertical edges which 
are not occupied by a dimer (magenta). The transfer matrix 
propagates the configurations in the positive y direction from 
row to row. Note the numbering of the “sites” (vertical edges): 
a fermion on site x may go to a: or a; -|- 1 in the subsequent 
row. A 3-plaquette is shaded. To enforce the presence of 
three dimers around this plaquette we need: (i) row j/ = 0: 
one fermion on edge 1; (ii) row y = 1'. one fermion on edge 1 
and one hole on edge 2 (or vice versa); (iii) row y = 2: one 
fermion on edge 2. The Eq. (18) is the associated expectation 
value. 

To treat the classical dimer model on the hexagonal 
lattice using a transfer matrix it is convenient to con¬ 
sider the “brick wall” version of the lattice, as depicted 
in Fig. 8. First, note that it is sufficient to consider the 
dimer occupations of the vertical edges - the information 
on all the other edges can be obtained using the hard¬ 
core constraints. The transfer matrix T then relates a 


y 



>1 


4 

‘ 

>3 
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>0 
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FIG. 9. Same as in Fig. 8, but now for a configuration with 
a 0-plaquette (shaded). To have this 0-plaquette, in terms of 
the fermionic representation, we need (i) row y = Q-. a. hole 
on edge 1 (il) row y = 1: two fermions on the edges 1 and 2, 
and (iii) row y = 2: a hole on edge 2. See Eq. (25). 

dimer configuration \ip) on one row y to the configura¬ 
tion on the row y +1 above. More precisely, T’lf/') is the 
linear superposition of all the configurations of the row 
{y -\- 1) which are compatible with [■!/)) at level y. The 
next step is to consider a single row of vertical edges, 
and to associate a (spinless) fermion Fock space to it: 
edges not occupied by a dimer carry the fermions, and 
edges with a dimer carry (fermionic) holes. In particular, 
the Pauli exclusion principle enforces the dimer hard-core 
constraint. We also note that the fermions (and their 
world lines) correspond exactly to the strings discussed 
in the main text. The x component of the flux density 
is simply related to density of vertical dimers, which is, 
in turn, simply related to the fermion density n (which 
is the same for each row): 

/ = 2 - 3n. (9) 

Thus, the particle number conservation of the transfer 
matrix enforces the flux conservation. If we note c\. the 
fermion creation operator on site x (see the numbering 
in Fig. 8), the transfer matrix can be shown to obey: 

^4 = (4 + cl+i)^. (10) 

TI vacuum) = [vacuum). (II) 

In other words, a fermion on site x should propagate to 
X OT {x + 1) in the line above. Performing a Fourier 
transform of Eq. 10 gives 

H = + ( 12 ) 

where is the Fourier transform of cl. This shows that 
T is a product of operators acting independently on each 
Fourier mode. The solution is [26]: 

f= n (l + e*Mcfc). (13) 

fcG [—'7r,7r[ 

From this one can find the commutation relations with 
annihilation operators: 

Ckf = fck (1 + e*4 , (14) 

Ca;f = f {C^ + Ca:-l) . (15) 
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In the following it will also be necessary to commute Cx 
(and cj.) and T in the reversed direction compared to 
Eqs. (10) and (15). The results are now infinite sums: 

4 ^’ = t{cI- 4+1 + 4+2 - 4+3 + •••). ( 16 ) 

~ i^x ^x—l “t” ^x —2 ^x—Z “t“ ' * * ) T, (ll") 

When the y dimension of the lattice goes to infinity, 
only the eigenvector of T with the largest eigenvalue 
in the given flux sector needs to be kept. The later is 
nothing but a Fermi sea |/) with Fermi momentum kp 
and density n = kpj-K. The corresponding eigenvalue, 
A{kp) = n-fcj,.<fe<fcj,.(l allows to compute the en¬ 

tropy per site, but its explicit expression is not needed 
here. 


Density of 3-plaquettes 

We start by the computation of psif), the density of 3- 
plaquettes. A corresponding hexagon is shaded in Fig. 8, 
and it is characterized by one fermion in x = 1 on the 
lowest row (thus associated with the projector cjci), one 


fermion in x = 1 and one hole in x = 2 in the second 
row (—>■ c|ciC 24 )j and, finally, a fermion at x = 2 in the 
third row (—>■ C 2 C 2 ). The density pa is thus 


r,{f\clc2Tc\ciC2clfc\ci\f) 
P3 — ^^- 


(18) 


(the factor 2 is due to the fact that there are two 
ways to put three dimers around an hexagon). The 
next step amounts to eliminate T by using the relations 
Fqs. (10),(15),(16) and (17). The result is 


P3 = 2(/|IiJ(c 2 -k ci)cic|4c2(c| c|)5i|/), (19) 


where we have defined: 

OO 

Dl = ^(-1)^4+., (20) 

x—0 

0 

Sr= Y. (-l)"c.+r. ( 21 ) 

X — — 00 

The correlator of Fq. (19) can be obtained, using Wick’s 
theorem, as the determinant of a 4 x 4 matrix M : 


Ma 


/ {D\{c2+Ci)) 

-((C2 +Ci)4) 

-((C2 +Ci) 4 ) 

V -((C2 -kCi)(cl +c\)) 


{blc) ( d | c 2 ) {bis,) \ 

-(ci4) (4^2) (44) 

-(ci4) (4^2) (44) 

-(ci(4 + 4)) -(c2(4 + 4)) ((4 + 4)4) / 


( 22 ) 


The two-point functions appearing above can be 
expressed using the correlator of the Fermi sea: 
Gx-y = (4cy) = and (4c^) = n. 

The correlations {blcy) or {cxS^) contain some infi¬ 
nite sums which can be evaluated using the sum rule: 

= n/2. The one appearing in Ma are 

{bid) = {blc 2 ) = {clSi) = (44) = n/2. The last 

one, {blSi), contains two sums which can also be per¬ 
formed exactly, leading to (i)J4) = 27 r[i+ 4 j(i^)] ' 
matrix Ma therefore takes the explicit form: 


Ma: 

Pa = 2det(Ma), 

([2-1- cos (utt)] n'^ — 2n + l) sin (utt) 

TT [cos (utt) -I- 1] 

(24) 

TT'^ 

Density of 0-plaquettes 


Ma 


( n nl2 n/2 
A n — 1 Gi 
A Gi n 
V 2A A A 


sin(n7r) 

27r[l+cos(n7r)] 

n/2 

n/2 

n 


\ 


(23) 


The density of 0-plaquette can be obtained in a similar 
way. The starting point is the following correlator (see 
Fig. 9): 


{f\c2clfcldclc2fdcl\f) 

{f\T^\f) 


(25) 


After commutating one T to the right and the other to 
the left we get: 

where G, = have set A = Gi — 1 -I- n. The 

quantity pa is finally obtained from the determinant of po = 2(/|(c2 -I- d)blc\dclc 2 Si{c\ + 4 ) 1 /)- (26) 
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As for Pa, we construct a matrix from the two-point con¬ 
tractions and the result is: 


/"-l "/2 "/2 slSSSfoll 

A n Gi nl2 

A Gi n n/2 

\ 2A A A n-l ) 


(27) 


Finally po is obtained by taking the determinant: 
po = det(Mo), 

cos (riTr) (cos (utt) -|- 1 ) -|- n^-n^ (n — 1 ) — 2 


—n sin (utt) 


cos (nTr) (n — 2) -I- 2n — 3 
TT (cos (nTr) -I- 1) 


-^ sin (nTr) (cos (nrr) — 1). 


(28) 


Note that we also computed pi using the same method 
and we checked that the sum rule pa = pi -I- 2po is satis¬ 
fied. 


Appendix H: Some variational states based on 
simple string arrangements 

In the main text we mentioned that some density pat¬ 
terns observed in the QMC simulations can be described 
by combining some particular “building blocks”, called 
S'-, H- and F-strings. Below we elaborate on this idea. 

An S-string is a static zigzag configuration, as depicted 
in Fig. 2 (main text). It corresponds to a row of 3- 
plaquettes. As an example, the classical star crystal can 
be viewed as a periodic arrangement of such S-strings 
at an average distance of 1.5 (measured in units of the 
distance between two nearest hexagon centers), noted 
therefore as a “classical Si .5 crystal”. i7-strings are set 
of configurations where, starting from a zigzag configu¬ 
ration, the string can move up by one on every second 
column. As for F-strings, they can move up by one on 
every column. 

Let us first suppose that we start with an isolated S- 
string, surrounded only by plaquettes carrying less that 
3 dimers. Switching on the kinetic term of the Hamilto¬ 
nian, a 3-plaquette located along the S-string is allowed 
to flip, with the constraint that none of its neighboring 
plaquettes has already flipped (the condition for still be¬ 
ing a 3-plaquette). We dub this constrained quantum 
system “F-string”. Note that, upon iterated flips, new 
3-plaquettes can be generated which were not bounded 
by the initial F-string, increasing the lateral extension of 
the set of 3-plaquettes. Such string configurations will 
play a role in high flux sectors, which have low string 
densities, but we do not consider further these extended 
chains in this present description. 

These constrained F-strings are interesting objects for 
themselves (see below). But their main interest here 


comes while considering their regular arrangements a 
distance d: indeed, for d < 3, correlated flips on two 
neighboring such strings may create 0 -plaquettes in be¬ 
tween, leading to an energy cost in the Ug > 0 part of 
the phase diagram. The QMC simulations show regions 
where the ground state dimer densities display such lin¬ 
ear zigzag arrays of 3-plaquettes, regularly spaced at a 
distance d, which we therefore call Sd and Fd crystals. 
Now, other nearby regions in the phase diagram (again 
in the “fan” region) show different patterns, such that 
only one over two of the 3-plaquettes is found to flip sig¬ 
nificantly. We call these configurations iJ-strings {H for 
“half’), and their associated regular arrangements the 
Hd crystals. Note that, being second neighbours in the 
zigzag chain, these plaquettes are free to flip, and we 
can already consider these id-strings as having a sim¬ 
ple resonating nature. An interesting feature is that 
the condition for not generating interstring 0 -plaquettes 
is now that d > d^™H — 5/2, instead of the above 
d > dp^^p = 3 for the F-strings. We therefore face an 
interesting competition between Fd and Fid crystals: the 
F-string has potentially a lower kinetic energy, but their 
dynamics can generate 0 -plaquettes at a shorter inter¬ 
string distance (as compared to the Fid crystals). 

When Vo gets large (and positive) we find that the 
ground state energy remains negative (even for U 3 > 0 ), 
meaning that the system can simultaneously gain some 
kinetic energy through resonances on 3-plaquettes while 
keeping a vanishing pg (in the limit Vg —>■ 00). Fig. 10 
shows some of the 0 -plaquette free states which can be 
obtained by stacking F-, H- and F-strings. For these flux 
sectors, the density patterns observed in the simulations 
approximately match those of these simple ansatze. 

We now discuss the energies of the Fd and Hd quantum 
crystals. 

According to whether d is integer or half integer, neigh¬ 
boring strings have parallel or anti-parallel zigzag config¬ 
urations. It is also easy to relate this distance d to the 
flux sector: in unit of the first distance between plaque- 
tte centers, one finds d = 3/(2 — /) (see Eq. 9 and note 
that n = 1/d). 

Let us, again, first consider an isolated iJ-string with 
independent resonant plaquettes. One then forms the 
tensor product of the individual chain ground states, and 
check whether this state leads to a useful variational ap¬ 
proximation. A necessary condition is that this state 
should not contain any 0-plaquette. Such plaquettes 
would appear in between two such chains if d is too short. 
As said above, this will not occur whenever d > 5/2 
(equivalently / > 0 . 8 ), giving an upper bound for the 
ground state energy. In the U 3 = 0 case, this leads to 
F(/) < —(2 — /)/ 6 , for / S [0.8,2]. The QMC simula¬ 
tions in the / = 0.8 sector gives an energy ^ — 0 . 22 , in 
rough agreement with the approximate — 0.2 value found 
for Ho/ 2 . 

Let us now analyze an isolated F-string. Start with 
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(f) SsHs crystal (g) HnFn crystal 

3 3 XT 



FIG. 10. Some variational states constructed by stacking S-, H-, and F-strings. (a) The staggered state which has flux 
density / = 2 and corresponds to zero strings, (b) The ideal plaquette state with / = 0, where all plaquettes of one of the 
three triangular sublattices are in the resonating state !■{_}•) + |\^). This corresponds to /f-strings at distance 3/2. (c)-(g) 
Different variational states built as tensor products of S-, H-, and T-string states in order to approximate certain crystals in 
the fan region for large vo- In these examples, the strings are at their minimal relative distances, such that the states are free 
of 0-plaquettes. The plaquette shading corresponds to the average number of dimers per hexagon, using the same color scale as 
in Fig. 3 of the main text (and in Fig. 5). In particular, 3-plaquettes are colored in dark blue, 2-plaquettes in light gray-blue, 
and 1-plaquettes in yellow. 


a configuration where each plaquette carries 3 dimers, 
and flip one of these. The flipped plaquette is still a 3- 
plaquette, and can therefore be later flipped back. But 
its two neighboring plaquettes now only carry 2 dimers, 
and are “frozen”. Repeating the flips on another 3- 
plaquettes allows one to span the full Hilbert space for 
such an isolated chain, containing constrained configu¬ 
rations of flippable or frozen plaquettes. For an open 
chain of length L the associated Hilbert space dimension 
is the Fibonacci number F/,. For a closed periodic chain 
the dimension is + Fl_i. It is possible to show an 
exact correspondence with the Hilbert space of so-called 
Fibonacci anyonic chain [38]. We numerically studied 
the quantum dimer Hamiltonian on a single chain and 
obtained for ua = 0 a ground state energy per plaque¬ 
tte ~ E_p = —0.6035605(9). The next step consists in 
building a tensor product of these states, which avoid 
0-plaquettes. This is possible if d > = 3, which 

means / > 1, giving upper bound E{f) < (2 — f)Ep/3. 
For / = 1 (Fa crystal), this gives an upper bound Ep/S 
which is slightly lower than that obtained above for the 
H ^/2 state. 

Back to the numerical results, we found that, for ua = 
0 and 2.4 < vq < 12, the ground state belongs to the 
/ = 4/5 sector, with a symmetry well described by that 
of the H 2.5 crystal. Note also that the above variational 
argument suggests that, at large uq, the sector / = 10/11 
would be be close in energy to that of / = 4/5. The 
corresponding state is an alternation of F- and id-strings 
at average distance d = 2.75 (see Fig. 10). We indeed find 
numerically that the energy in the sector / = 10/11 falls 
below that of / = 4/5 when vq > 200, and displays the 
expected El — F density pattern. 


More complex states appear to be selected for vq > 
12. Some of these states can be qualitatively understood 
by introducing correlations among the strings. A simple 
family of building blocks consists in forming n coupled F- 
strings. There, plaquette flips are not only constrained 
by the state of the two neighbouring plaquettes along the 
string, but also along the perpendicular direction. One 
then form a tensor product of these blocks of correlated 
strings to construct a state of the whole system. This 
leads to the infinite series of fluxes f{n) = (4n-|-2)/(5n-|- 
1 ), running from / = 1 (single strings, F 3 crystal) to / = 
4/5 (all strings coupled, the H 2.5 case). We numerically 
found that the flux f = 6j7 compatible with 4 coupled 
strings indeed replaces the / = 4/5 sector as the ground 
state for uq > 12 . 

Appendix I: Transition from the / = 1/2 sector to 
the fan region 

In the main text, a variational result for the transi¬ 
tion line between the S' 2 -crystal (/ = 1/2) and the fan 
region (1/2 < / < 2) was stated for large vq- We dis¬ 
cuss it here in a bit more detail. As an ansatz for the 
S 2 ground state, we employ |... SSSSSSSSSSSS ...), 
a regular arrangement of static F-strings with distance 
ds-s = 2 between neighboring strings. It has an en¬ 
ergy of E/N = V 3/2 per hexagon. We compare it with 
a state that has an infinitesimally increased flux den¬ 
sity (one string less such that AFy = 3). For the cor¬ 
responding ansatz state, some F-strings are replaced by 
F-strings such that we change the configuration above 
to I... SSHSHSHSHSS ...) with (average) string dis¬ 
tances ds-s = 2 and ds-H = 2.25. To keep the total 
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system size constant size, note that we need to replace 
five ^'-strings by four iJ-strings. The resulting energy 
difference per column of hexagons is 

AE/L, = -5v3 + 4(^-^-+3'^y 

which vanishes at vs/t = —1. This shows that a simple 
S 2 crystal gets destabilized with respect to i7-string in¬ 


sertion for vs/t > —1. This simple variational argument 
thus predicts that the transition to the fan region occurs 
(for vq —> 00 ) at va/t « —1, indeed not far from the 
observed value. 
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